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13.  ADM  RAC  T 


An  error  analysis  is  presented  for  Gaussian  elimination 
when  the  matrix  is  arbitrarily  sparse.  Error  analyses  for 
elimination  on  band  matrices  and  full  matrices  follow  as 
special  cases. 
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1.  INTRODUCTION 


Since  the  direct  solution  of  systems  of  linear  equations 
by  elimination  is  now  well  understood  when  the  matrix  is  full 
or  band,  attention  has  turned  to  the  problem  of  elimination 
of  systems  when  the  matrix  is  arbitrarily  sparse.  When  the 
sparse  matrix  arises  from  the  discretization  of  an  elliptic 
partial  differential  equation,  it  has  a  special  structure  and 
special  properties  which  make  iterative  methods  attractive. 
Even  here,  however,  one  can  often  find  an  ordering  of  the 
equations  which  makes  elimination  competitive  (George [3]) . 

But  all  too  frequently  the  matrix  has  an  arbitrarily  sparse 
structure,  e.g.  in  network  analysis  and  econometric  problems. 

Here  we  present  an  error  analysis  of  (point)  Gaussian 
elimination  when  the  matrix  is  arbitrarily  sparse.  The  anal¬ 
ysis  is  presented  in  terms  of  the  structure  of  the  elimina¬ 
tion  graphs  arising  in  the  elimination  process.  The  only 
assumption  we  make  is  that  the  leading  principal  minors  are 
non-zero,  i.e.  that  the  (1,1)  entry  in  each  reduced  matrix 
is  non-zero  (cf .  Forsythe  and  Mo.ler  [2]  ,  pp.  27-36)  . 

Since  the  graphs  of  band  and  full  matrices  are  special 
elimination  graphs,  the  error  analyses  for  band  and  full 
matrices  follow  as  special  cases. 


2.  MATRICES  AND  GRAPHS. 

Let  M  be  an  n*n  matrix  with  non- zero  leading  principal 
minors;  hence  M  is  non-singular  and  the  LU  decomposition  of 
M  exists  and  is  unique.  If  M  is  symmetric,  then  M=LDL  . 

Rose  (4]  associated  an  undirected  graph  with  a  sym¬ 
metric  matrix  and  interpreted  the  LDLt  factorization  graph- 
theoretically  by  undirected  elimination  graphs.  Bunch  and 
Rose  [1]  made  the  extension  of  the  association  of  directed 
graphs  with  general  square  matrices  and  of  the  interpreta¬ 
tion  of  the  LU  factorization  by  directed  elimination  graphs. 

The  directed  graph  of  M,  G(M ,  with  ver Lex  set 
X  and  arc  set«,is  defined  as  follows;  a  vertex  x^eX  is 
associated  with  row  i  of  M,  and  (x^,Xj)e«  (an  arc  from  x^ 
to  x.  is  in  G)  if  and  only  if  m.  .j^O  and  i^j.  The  vertices 

J  J 

X  are  regarded  as  ordered;  i.e., 

[a  rt 

Let  the  matrix  M  be  written  as  M= 

c  B 

is  lxl,  r  and  c  are  (n-1) xl  and  B  is  (n-l)x(n-l).  Then 
the  first  step  of  the  LU  factorization  of  M  can  be  written 


If  G(M)  is  the  directed  graph  of  M,  the  elimination 


graph  is  obtained  from  G  by  deleting  y  together  with 

its  incident  arcs  and  adding  an  arc  (x,z)  whenever  there 

exists  a  (directed)  x,  z  path  of  length  2  containing  y. 

Gy  is  the  graph  of  the  matrix  obtained  by  eliminating" 

the  variable  corresponding  to  y  in  Gaussian  elimination? 

e.g.,  G  is  the  graph  of  B-crVa.  The  accidental  ere- 
X1 

ation  of  zeros  during  the  elimination  process  is  ignored. 
For  a  more  detailed  discussion,  see  Bunch  and  Rose  [ 1] , 
Section  2. 


Let  G^ , . . , »Gn_^  be  the  sequence  of  elimination 
graphs  defined  recursively  by  GQ=G(M)  and  G^-(G^_^)x^. 
Let  be  the  number  of  elements  in  the  set  Let 
ri=l (VcXi.1:!xi,y) e«1_1,G._1=<x._1,oi_1) )|  and 

C^l  {yeXi-;L:  (y,xi)e«ri-1,Gi-1=(Xi-1,^i-;i)  }|  be  the  out- 

degree  and  in-degree,  respectively,  of  vertex  x^  in  the 
elimination  graph  G^_^. 

Note  that  r^+1,  c^+1  is  the  number  of  non-zero  ele¬ 
ments  in  the  first  row,  column  of  the  reduced  matrix  of 
order  n-i+1,  i.e.  the  reduced  matrix  whose  graph  is  G^_j, 
Define  e^^l  if  there  is  an  arc  from  x^  to  x^  and 
eki”l  there  is  an  arc  from  x^  to  x^  in  Other¬ 

wise,  let  e^=0  and  e^^=0.  We  count  a  division  as  a  mul* 
triplication  and  a  subtraction  as  an  addition. 


4 


When  M  is  symmetric,  then  r.=c. =d.  ,  and  Rose  [4]  has 

t  n-1 

shown  that  the  factorization  M=LDLr  requires  E  d. (d.+3)/2 

i=i  1  1 

n-1  1 1 

multiplications  and  E  d.  (d.+j.)/2  additions,  while  the 

i=1  1  1  n  , 

t  n  1 

backso.lving  of  I.DL  x=b  requires  n+2  E  d.  multiplications 

i--l  1 

n-1  1  1 

and  2  E  d.  additions. 
i”l  1 

When  M  is  general.  Bunch  and  Rose  [1]  have  shown 
that  the  M«LU  factorization  requires 


n-1 

E  (r.+l)c.  multiplications  and 
i=l  1  1 


n-1 

E  r.c.  additions,  while  the 
i“l  1  1 

backsolving  LUx=b  requires 
n-1 

n+  E  (r.4c. )  multiplications  and 
i~l  1  1 

n-1 

E  (r.+c.)  additions. 
i=i  1  1 


Examples :  Let  M  be  an  n*n  matrix. 

(1)  If  M  is  symmetric  and  full,  then  d^=n-i  and  solving 
Mx--LDLtx=b  requires  n+  E  d^(U^+7)/2=  |n3  +  |n2  -  |n  multiplica¬ 
nd  1327 

tions  and  E  d.  (d.+5)/2  =  ?n  +  n  -  zn  additions. 
i~X  ^  0  v 


(2)  If  M  is  symmetric  and  band,  with  bandwidth  2m+l, 


then  d^=m  for  l£i£n-m  and  d^=n-i  for  n-m+l£ij<n-l.  Solving 

i.  127  1325 

Mx=LDL  x=b  requires  (|m  +  ^m  +  l)n  -  ^m  -  2m  -  -m  mul- 

12  5  13  32  7 

tiplications  and  (^m  +  ^m)n  -  additions. 

(3)  If  M  is  general  and  full,  then  r.=n-i=c.  and 

13  2  1 

solving  Mx=LUx=b  requires  n+  E  (r . c.+2c.+r. )=  -^n  +  n  -  ^n 

n— 1  13  12  5 

multiplications  and  Z  (r^c^+c^+r^)  =  ^n  +  in  -  ^n  additions. 


(4)  If  M  is  general  and  band,  with  bandwidth  2m+l,  then 
r^=m=c^  for  l£i£n-m  and  r^=n-i=c^  for  n-m+l£i<n-l  when  no 
pivoting  is  needed,  while  c^=m  for  l^i^n-m,  c^=n-i  for 
n-m+l£i<n-l,  r^^^m  for  l<i<n-2m,  r^n-i  for  n-2m+l<i_<n-l , 

when  partial  pivoting  is  used.  Solving  Mx=Ll3x=b  requires 

2  2  3  2  4 

(m  +3m!l)n  -  ^m  -  2m  -  ^ m  multiplications  and 

7  7  3  3  7  5 

(m  +2m)n  -  ^m  -  ^m  -  additions  when  no  pivoting  is 
2  13  3  2  11 

needed,  and  £(2m  +4m+l)n  -  g— m  -  4m  -  g— m  multiplications 
7  133  72  4 

and  <(2in  +3m)n  -  ^m  -  jin  -  |m  additions  when  partial 


pivoting  is  used. 
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3.  ERROR  ANALYSIS  OF  TRIANGULAR  FACTORIZATION. 


n 


Let 


n 


=  max  |x.j  ,  { |  M{  H  max  Z|xn..|, 
l<i<n  1  l<i<n  j=l  13 

..  n  . 

1 1  x|  j  .  =  E  |  x.  |  ,  |  |  Mj  |  ,  =  max  l  |  m.  .|  .  If  M=M  ,  then 
1  i-1  1  1  l<j<n  i=l  13 

|  |  M|  |  x-|  {  M|  |  eo. 

Due  to  the  finite  precision  arithmetic,  we  obtain  the 

triangular  factorization  of  a  perturbation  of  M,  i.e.  LU=M+F. 
n  \ 

Let  M  =  M.  Then  the  elimination  is  defined  sequential¬ 
ly  for  lj<k£n-l  by  the  following:  given  f  let  L^  be  de¬ 
fined  by  i,|k^=l  for  ^ik^ ^mik^ //rakk^  f°r  n' 

otherwise,  then  M(k+1 }  ~f  l  ( (L(k) ) -1M(k) )  ,  i.e.  mfk+1)  = 


r 


n  00 

mii 

V.  13 


for  k+l<i<n,  j=k 
for  k+l<i ,  j  <n 
otherwise. 


Then  U=M*n)  and  L=L(1)L*2* * • •L(n_1) ,  i.e. 


LiJ 


0 

1 

L<3> 


^i'-i. 


for  i< j 
for  i=j 
for  i>j 


Following  the  notation  of  Forsythe  and  Moler  [2],  §21, 
let  u  be  the  unit  round-off,  e.g.,  u=  for  rounded 

operations  in  t  digit,  base  8  arithmetic.  Then 


*TOik)/mkk^  (1+5ikeik*  '  where  !  5i)cl  — u '  or 


0  =  lnil^)'£ik)lnkk)+eik>  '  Where  e4V>Sm<V>64Ve4V  for 


ik  ik  °ikcik 


k+l<i£n ,  and  for  k+l<i,  -j<n: 


-S*1’  ■  (I»ij,-4,'”)‘j)<11{i3cik0k:»/<l+5i:eikeicj>' 

where  <$^-,=0,  I  5ii.!  Jl11  for  2i*-5n'  iSijilu  for  j<n,  or 


Jk+1)  _  (k)  g  (k)  (k)  (k) 

mij  =  mij  *ik  mkj  +Eij  ■ 


where  e<: f  f  «ijeikekj‘mij+:1  ’  6ijeik°kj  ‘  Let  eij)=0 


otherwise . 

Let  |  |  <t  for  all  i,  k  and  |  |  <cj  for  all  i,j,k. 

Then  I  ei^  I  laueik  for  k+l^ifn  and  |  I  <.(t+l)  aueikekj  for 
k+l£i ,  j£n. 

Let  F^  =M*k+1*-M^+L^M^  for  l<k<n-l; 


•  <* 
l-l 

n-l  ...  u  (t+1)  o  E  e..e.  . 

then  F=  E  FW  and  |F..|  </  k=l  1  3 

k=l  ^  j-1 

“(T+1)\^eikekj+uc,eij 


for  i£j 


for  i>j 
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(1-1  ji-1  j-1  n  i-1 

Thus,  i|F||^ua  max)  E  e..+  (t+l)  E  E  e..e..+  E  E  e..e,  . 

lii<n\j=1  13  j=l  k=l  1K  K]  j-i  k=l  1K  *3 


r  n  3  1-1  n  3-1 

and  |  J  F|  j  ,  <ua  max  1  E  e..+  (T+l)  E  E  e.,e,  .+  E  E  e.,e 

1<j<nU=j+i  3  i=2  k~l  1K  *3  i=j+l  k=l  1K 

If  M=Mt  and  LDLt=M+F ,  then  F=Ft ,  for  all  i,j, 


i-1  j-1 

|Fu|<u(t+1)o  £  sik,  |Fij|<u!T+l)o  Z  eike.k+uoei:j  for  i>j, 
K—±  K— J. 


±  —  ± 

I  Fijl  £u(t+1)  0  E  eijcejk+uaei j  for  i< j  ,  and  ||  F|  |  ro=|  |  p|  (  L= 

}c — 1 

n  ( n  i-1  j-1  i-1 

max  E  |F..|<ua  max  I  E  e..+  (T+l)  E  E  e.,e..,+  E  e.,+ 
l<i<n  j=l  13  "  l<i<n\ j=l  30  j=l  k=l  1K  jK  k=l  1K 

L 


n  i-1 


2  **  eikeik  \  ' 

j=i+l  k=l  lk  3K  / 


* 
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4.  ERROR  ANALYSIS  OF  BACKSOLVING. 


Solving  Ly=b,  where  L  is  unit  lower  triangular,  we  ob¬ 
tain  the  exact  solution  of  (L+SL)y=b.  Here  y.j=b^  and 

yi-fi(-iily1-ii2y2 - ^,1-^-1^)  = 

(“s'il(l+5il)yi“* '  ‘“ai'i-l^1+5i,i-l^yi-l+bi^1+<Sii*  '  where 

i-1 

6n=0,  |6^j£u  for  2£i_<n,  |  5I. Olu (1+  Z  eki^eii  f°r  2£i£n, 

k~2 

i-1 

I  6 .  .1  <1.01u(l+  Z  e..)e.  .  for  2<j<i<n.  Then 

13  ”  k=j+l  *3  13  ~  ~ 

n 

I I  6L|  I  ,  <1.01u  maxi  l .  .1  max  {  Z  |6..|}< 

i,j  13  l<j<n  i==j  13  ” 

1.01uTmax{icn {c,+l) ,  max  l^c . (c .+1) +1] }  since 
1  1  1  2<j<n-l  1  3  3 

n  n  i-1 

Z  I  <5.  J  <1.01u  Z  (1+  Z  ev,)e.  »=1.01uc,  (c,+l)/2  and 
i=l  11  “  i=2  k=2  Ki  11  1  i 


Z  iSi^fl.Olu  {1+  Z  (1+  Z  ^ek;.)eij}  ^l.OluI^  (Cj+l)+ll 


i=j+l  k=j+l 


1  3 


for  2<j<n-l,  while  1  I 5L|  I  <1.01ut  max  {e. ,+  Z  e.,e., 

l<i<n  11  k=2  kl  11 

i-1  i-1  n-1  n-1 

+  Z  (1+  Z  e.  .)e.  .+1} <1 . 01ux{l+c, -1+  Z  c .+1}=1. OIut (1+  Z  c.) 
j=2  k=j+l  k3  13  j=2  3  j=l  3 


Solving  Ux=y,  where  U  is  upper  triangular,  we  obtain 
the  exact  solution  of  (U+6U)x=y.  Here  xn"f  £(yn/unn)  = 


V1"nn(UUl  ana  xi=£* 


i,i+l  i+1 


xi+i — uinVn 
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“ui,i+l(1+gi,i+l)xi+l  *'"uin  (1+gin)xn+^i 


for  2<i<n, 


where  $^n=0,  J6!.|<u  for  l<i<n-l,  !  5 ^ 4 1 for  l<i<n. 


^•5nll1*01u^1+  s  evr.^e-ir>  -or  l<i<n“l,  and 
1  k=i+l  Kn  in  - 


1 j|  <1-01x1(2+  £  ev.)e..  for  l<i<j<n-l.  Then 

13  k=i+l  *3  13  -  - 


I  I  6UI  I  l^axlu..!  max  {|6J.|  +  £  |6..|}  < 
1  i,j  13  l<j<n  33  i=l  13  ~ 


j-1  j-1  n-1  n-1 

1.01uamax{  max  [2+  £  (2+  £  e.,)e..],  1+  £  (1+  £  e.  )e.  }, 

1< j<n-l  i=l  k=i+l  K3  13  i=l  k=i+l  xn  in 


while  1 1  6U|  l^maxlu.j  max  {l6!.j+  £  |  <$..|  }  < 
i,j  J  l£i<n  Xl  j=i  13 

n-1  j-1  n-1 

l.Oluo  max  {2+  £  (2+  £  e..)e.  .+  (1+  £  e.  )e.  }  < 

l<_i<_n-l  j=i+l  k=i+l  *3  13  k=i+l  kn  in 


l.Oluo  max  {2+2r. -2e.  +^r. (r. -l)+e.  }  < 
l<i<n-l  1  m  2  1  x  in  - 


l.Olua  max  {2+£r. (r .+3) } . 
l<i<n-l  1  1  1 


When  U=DLfc,  ri=ci"^i  an<^  we  obtain  (D+6D)z-y  and 

(Lt+ALt)x=z.  Here  z^=fS,(yi/d^i)=yi/[di^(l+6^i)  ]  where  I  5 ill  ^.u 
for  l£i<n,  so  |  )  6D|  |  3=|  ]  6d)  )  ^<u  max|  d^J  =up,  where  p=max|  <3^i|  . 
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1 1  ALfc|  I  ,<1.01ut  max{  max 
11  i-  l<j<n-i 


j-1  j-1 

[1+  E  (1+  E  e. 
i=l  k=i+l  K3 


n-1 

E 

i=l 


n-1 
E  e 
k=i+l 


knein 


} 


and  1 1  1^1.0  lux 


max  {1+id. (d.+l) } . 
l<i<n-l  1  1 
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5.  ERROR  ANALYSIS  OF  ELIMINATION . 


Solving  Mx=b  in  finite  precision,  we  obtain  the  exact 

solution  of  (M+6M) x=b,  where  6M=F+ (6L) U+L(5U)+ (6L) (SU) ,  or 

when  M-Mfc,  5M~F+LD  { ALt)  +  [L  { 6D)  +  ( <SL)  D+  ( <$L)  (6D)]  [L^+AL^J  . 

i-1 

||l[|.<1+t  max  c.,  |  J  L|  |  <1+t  max  E  e .  .<1+ (n-1)  t  , 

3  2<i<n  j=l  13~ 

n 

1 1  011 -  <0(1+  max  E  e.  .)<na,  [  |  U  [  |  a  ( 1+  max  r.), 

j=i+1  3  ”  l<j<n-l  3 

l|Lt||1“l|L||„,  II^IL-llLHi,  and  IlDll^llDll^p. 


(t+1) 


{it 

E  e.  .+ 
i=j+l  13 

f  j  i-1  n  j-1 

li=2 


+  T 


1  11 
max  (1+pC . (c  .+1) ] [1+  max  E  e..] 

l<j<n-l  3  3  l<i<n-l  j=i+l  13 


j-1  j-1 

(2+  E  ek.)e..], 
k=i+l  *3  13 


(j  * 
max  [2+  E  1 
!<:  j<n-l  i=l 

n-1  n-1  *\ 

1+  E  (1+  E  e.  ) e.  }  (  1+t  max  c. 
i=l  k=i+l  Kn  in;  v  l<j<n-l  3 

+1.0 Iut  max  [1+ic. (c.+l) ]}  >  .  j  |  5M|  |  can  be  bounded 
1< j<n-l  1  3  3 


similarly.  Similar  bounds  are  obtainable  when  M=M^. 
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6.  THE  BAND  AND  FULL  CASES. 

When  M  is  a  band  matrix  with  bandwidth  2m+l,  then  e^j=0 
for  |i-j|>m;  under  Gaussian  elimination  with  partial  pivoting, 
e^=0  for  i>j+m  and  for  i<j+2m. 

With  no  pivoting,  |  |  F|  |  co<au[m+ (T+l)m'  ] ,  1 1  L|  |  ^1+mx, 

1 1  6LI  I  ^1.01ux(|m2+im+l] ,  |  1  U|  |  ^ (m+l)a,  and  |  |  5U|  |  oa< 
1.01ua[^m2+|m+2]  ;  thus  [  |  6M|  |  ^1.01uc{xm3+|(5x+3)m2+ 

|(7x+S) +T+2+1 .  Olux (^m4+m3+|m2+|m+2) } . 

If  M  is  diagonally  dominant,  then  t<1,  a<2,  and  |  |  6M|  | 
1.01u{2m3+  8m2+12m+4+l. Olu (^m4+2m3+^-m2+5m+4)  }  . 

With  partial  pivoting,  x«l,  cK22m-3'- (rn-1)  2m*"2  (Wilkinson [5] )  , 
I  hi  |  00<uc[5m2+4m]  ,  |  |  L|  |  ^l+mx  ,|  |  6L|  |  ^1.01u(im2+|m+l]  ,  |  |  u|  |  m  < 
(2m+l)o  ,  and  |  j  6u|  |  ^l.Olua  [2m2+3m+2]  ;  thus  |  |  6M|  |  to  _< 
1.01ua{3m3+-^m2+— |m+3+1.01u(m4+|m3+|m2+4m+2)  } . 

When  M  is  full,  e.  .=1  for  all  i,j.  Thus  | | f| |  < 

13  00  _ 

ucr{n-l+  (t+1)  (|n2-“n)  }=iuon‘[  (t+1)  n-x+1] ,  |  [  I.|  |  ^<1+  (n-1)  t  , 

lluH^no,  1  1  6L|  I  ^<1.0 lux  [n  (n-1)  /2+1]  ,  |  |  6U|  |  oo<1.01uon(n+l)/2, 

and  |  |  6M|  |  „< 1.0 lug  [xn3+n2+  (2-t)  n+  1  .Olux  (n4+2n2+2n)  /4]  . 

With  partial  or  complete  pivoting,  x=l  and  |  |  6mJ  [  m  £ 
1.01ug[n3+n2+n+l .01u(n4+2n2+2n) /4] . 
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7.  REMARKS. 

From  §§3-5  we  see  that  the  error  matrix  arising  in  sparse 
elimination  depends  on  the  fill-in  occurring  during  elimina¬ 
tion.  In  the  band  case  we  know  a  priori  the  structure  of  the 
matrix.  However,  in  the  arbitrarily  sparse  case  we  must  view 
elimination  in  the  context  of  optimal  ordering,  i.e.  an  or¬ 
dering  of  the  equations  so  that  fill-in  is  minimized.  Given 
a  fixed  ordering,  our  analysis  here  bounds  the  errors  occur¬ 
ring  in  the  computation. 
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